#!/usr/bin/perl -w
use strict;
use Bio::SeqIO;
use Bio::Seq;
use Bio::PrimarySeq;
use Bio::Tools::Run::StandAloneBlast;
#Copyright 2011. Jason Weirather.
#This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License v3 or any later version.
my $mcloutfile = shift @ARGV;
my $indexfile = shift @ARGV;
my $outfile = shift @ARGV;
my $fastafile = shift @ARGV; #we need the sequences of course.

my $minimum_group_size = 8;  # be sure if you change this to also change it in the automated database annotation program cluster_conserved_sequences.pl

my $ival = shift @ARGV;
#read in the base results
open(IIN,"$indexfile") or die "couldn't open that file\n";
chomp(my @indexlines = <IIN>);
close IIN;
my %map;
my $firstline = shift @indexlines;
my $secondline = shift @indexlines;
if(!($firstline=~/^#/) && !($secondline=~/^#/)) { die "format error\n"; }
foreach my $indexline (@indexlines) {
  my ($index,$name) = split(/\s+/,$indexline);
  $map{$index} = $name;
}

open(FASTAIN, "$fastafile") or die " could not open $fastafile the fasta file\n";
my $seqin = Bio::SeqIO->new(-fh => \*FASTAIN,
                             -format => 'fasta');
my %seqs;
while(my $seq = $seqin->next_seq()) {
  my $seqid = $seq->id;
  if($seqid=~/^(\S+)\s+/) { $seqid = $1; } #first nonwhitespace.  you may have to do this for the other first but i think mcl is stripping it too.
  $seqs{$seqid} = $seq->seq;
}
close FASTAIN;

open(MIN,"$mcloutfile") or die "couldn't open that file\n";
chomp(my @mlines = <MIN>);
close MIN;
my $conc = '';
while(my $line = shift(@mlines)) {
  $conc .= "$line ";
}
if($conc=~/\(mclmatrix\s+begin\s+(\S.*)\s+\)/) { $conc = $1; } else { die "couldn't read mcl output\n"; }
my @outlines = split(/\s*\$\s*/,$conc);
open(OBF,">$outfile") or die;
foreach my $line (@outlines) {
  my $classlabel = '';
  if($line =~/^(\d+)\s+(.*)$/) {
    $classlabel = $1;
    $line  = $2;
  } else { die "still couldnt read\n"; }
 my @indecies = split(/\s+/,$line);
 if(scalar(@indecies) >= $minimum_group_size) {
  open(OTF,">$outfile.$classlabel.fasta");
  ###this is where we try to figure out orientation and fix it.  TOUGH
  #foreach my $index (@indecies) {
  my $firstseq;
  my $isProtein = 0;
  my $longest = 0;
  my $longest_index = 0;
  my $count = 0;
  foreach my $index (@indecies) { 
    if( $seqs{$map{$index}} =~ /[RDQEHILKMFPSWYV]+/) { $isProtein++; } 
    if(length($seqs{$map{$index}}) > $longest) {
      $longest = length($seqs{$map{$index}});
      $longest_index = $count;
    }
    $count++;      
  }
  $firstseq = Bio::PrimarySeq->new(-seq => $seqs{$map{$indecies[$longest_index]}}, -id => 'fir'); 
  for(my $i = 0; $i < scalar(@indecies); $i++) {
    #print "$index  $classlabel\n";
    if(!exists($seqs{$map{$indecies[$i]}})) { die "what?\n"; }
    my $curseq = Bio::PrimarySeq->new(-seq => $seqs{$map{$indecies[$i]}}, -id => 'cur');
    if($i != $longest_index) { #compare to the first one.  The higher scoring orientation wins
      my @params;
      if($isProtein == 0) { @params = (program => 'blastn'); }
      if($isProtein > 0) { @params = (program => 'blastp'); }
      my $factory = new Bio::Tools::Run::StandAloneBlast(@params);
      my $aln = $factory->bl2seq($firstseq,$curseq);
      eval {
        my $strand = $aln->next_result->next_hit->next_hsp->strand('hit');
        if($strand == -1) { # fix it to be like the first strand
          $seqs{$map{$indecies[$i]}} = $curseq->revcom()->seq;
        }
      }
    }
    print OBF ">$map{$indecies[$i]}\t$classlabel\n$seqs{$map{$indecies[$i]}}\n";    
    print OTF ">$map{$indecies[$i]}\t$classlabel\n$seqs{$map{$indecies[$i]}}\n";    
  }
  close OTF;
 }
}
close OBF;
